//
// CUDA Sagan Simulator - a CUDA version fo the sagan Simulator Software
// version 1.3 - Alpha 01
// by TTYMan
//
// para compilar:  nvcc SS.cu -o sagan -I/home/ttyman/NVIDIA_GPU_Computing_SDK/C/common/inc -arch sm_21
// para plotar grafico das energias no gnuplot: plot 'energy.log' using 2 title 'ETotal', 'energy.log' using 3 title 'EK', 'energy.log' using 4 title 'EP'
//

#include <stdio.h>
#include <string.h>
#include <cutil.h>
#include <cuda.h>
#include <cuda_runtime.h>

/*
#if __CUDA_ARCH__ < 130
typedef float REAL;
#else
typedef double REAL;
#endif
*/

typedef double  REAL;
typedef double2 REAL2;
typedef double4 REAL_CUDA;

//cutil.h ==> /home/ttyman/NVIDIA_GPU_Computing_SDK/C/common/inc

#define BSIZE 256
#define BUFFER_MAX 2048      // buffer size for reading each line of initial conditions file.
#define offset      4        // this offset space is used to define particle position information (float4 CUDA datatype):
                             //   0 (x) - Pos X ==> array[index*offset+0]
                             //   1 (y) - Pos Y ==> array[index*offset+1]
                             //   2 (z) - Pos Z
                             //   3 (w) - Mass
                             // The velocity information is stored ina float3 CUDA datatype, as follows:
                             //   0 (x) - Vel X
                             //   1 (y) - Vel Y
                             //   2 (z) - Vel Z
                             //   3 (W) - Not used (reserved for future use).
                             // this approuch is useful to take advantage of CUDA float4 datatypes.


#define SHARED_KERNEL 96
#define SHARED_MEMORY_SIZE 4096

#if __CUDA_ARCH__ <300
#define BLOCKS_PER_MP 8
#else
#define BLOCKS_PER_MP 16
#endif

// declaracoes memoria device

// device memory for particles array
__device__ REAL_CUDA*  d_ParticleArrayPosition;
__device__ REAL_CUDA*  d_ParticleArrayVelocity;

// device memory for energy monitoring
__device__ REAL2* d_EnergyArray;
// d_EnergyArray.x = Potential Energy for the particle
// d_EnergyArray.y = Kinetic Energy for the particle


__device__ REAL_CUDA BodyBodyInteraction(REAL_CUDA ParticleIPosition, REAL_CUDA ParticleJPosition, REAL_CUDA Acceleration, REAL SofteningFactorSQR)
{
    REAL XDist;
    REAL YDist;
    REAL ZDist;
    REAL DistanceSQR;
    REAL Rinv;
    REAL Rinv3;

	// [1 FLOP]
	XDist = ParticleIPosition.x - ParticleJPosition.x;

	// [1 FLOP]
    YDist = ParticleIPosition.y - ParticleJPosition.y;

    // [1 FLOP]
	ZDist = ParticleIPosition.z - ParticleJPosition.z;

    // [5 FLOP]
	DistanceSQR = (XDist * XDist) + (YDist * YDist) + (ZDist * ZDist);

    // [4 FLOP]
    Rinv = 1.0 / sqrt( DistanceSQR + SofteningFactorSQR );

    // [2 FLOP]
    Rinv3 = (Rinv * Rinv * Rinv);

    // [9 FLOP]
    Acceleration.x -= ParticleJPosition.w * XDist * Rinv3;
	Acceleration.y -= ParticleJPosition.w * YDist * Rinv3;
	Acceleration.z -= ParticleJPosition.w * ZDist * Rinv3;

	// [5 FLOP]
	// The 'w' part is used to save the Potential Energy value calculated for this pair of particles.
	Acceleration.w -= (ParticleIPosition.w * ParticleJPosition.w) / sqrt( DistanceSQR + SofteningFactorSQR );

    return Acceleration;
}

__device__ void FillSharedMemoryBuffer(REAL_CUDA* d_ParticleArrayPosition, unsigned int N, unsigned long tile, unsigned long tilesAmount, int tileWidth)
{
	extern __shared__ REAL_CUDA sh_ParticleArrayPosition[];

	// last tile?
	if(tile == tilesAmount-1)
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			if(tile*tileWidth+i < N)
			{
				sh_ParticleArrayPosition[i] = d_ParticleArrayPosition[tile*tileWidth+i];
				//sh_ParticleArrayPosition[i] = d_ParticleArrayPosition[tile+i];
			} else break; // avoid out of bound access on shared memory buffer.
		}
	}
	else
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			sh_ParticleArrayPosition[i] = d_ParticleArrayPosition[tile*tileWidth+i];
		}
	}

}

__device__ REAL_CUDA TileProcessing(REAL_CUDA* d_ParticleArrayPosition, unsigned int N, REAL_CUDA Particle, REAL_CUDA Acceleration, REAL SofteningFactorSQR, unsigned long tile, unsigned long tilesAmount, int tileWidth)
{
	extern __shared__ REAL_CUDA sh_ParticleArrayPosition[];

	// last tile?
	if(tile == tilesAmount-1)
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			if(tile*tileWidth+i < N)
			{
				Acceleration = BodyBodyInteraction(Particle, sh_ParticleArrayPosition[i], Acceleration, SofteningFactorSQR);
			} else break; // avoid out of bound access on shared memory buffer.
		}
	}
	else
	{
		for(unsigned int i=0; i<tileWidth; i++)
		{
			Acceleration = BodyBodyInteraction(Particle, sh_ParticleArrayPosition[i], Acceleration, SofteningFactorSQR);
		}
	}

	return Acceleration;
}



__global__ void MainKernel(REAL_CUDA* d_ParticleArrayPosition,
                           REAL_CUDA* d_ParticleArrayVelocity,
                           REAL2* d_EnergyArray,
                           unsigned int N,
                           REAL SofteningFactor,
                           REAL TimeStep,
                           int tileWidth)
{
	// map shared memory buffer
	extern __shared__ REAL_CUDA sh_ParticleArrayPosition[];

	// global thread ID: unique ID for the particle.
	unsigned int gThID = blockIdx.x * blockDim.x + threadIdx.x;

	if (gThID>N) return;

	REAL_CUDA Particle = d_ParticleArrayPosition[gThID]; //global memory access!

	REAL_CUDA Acceleration = {0.0, 0.0, 0.0, 0.0};

	// save some FLOPs squaring epson factor here, avoiding to do this inside the tile processing loop.
	REAL SofteningFactorSQR = SofteningFactor * SofteningFactor;

	// how many tiles will be used?
	unsigned long tilesAmount = N / tileWidth + (N % tileWidth == 0 ? 0 : 1);

	for(unsigned long tile=0; tile < tilesAmount; tile++)
	{

		// the first block thread will be used to fill the shared memory buffer.
		if(threadIdx.x==0)
		{
			FillSharedMemoryBuffer(d_ParticleArrayPosition, N, tile, tilesAmount, tileWidth);

		}
		__syncthreads(); //make sure shared memory buffer was updated for all threads.

		Acceleration = TileProcessing(d_ParticleArrayPosition, N, Particle, Acceleration, SofteningFactorSQR, tile, tilesAmount, tileWidth);

		//__syncthreads();
	}

	//__syncthreads();


	//d_EnergyArray[gThID].x = 1.0;
	//d_EnergyArray[gThID].y = 1.0;

	// updating the velocities
	d_ParticleArrayVelocity[gThID].x   += Acceleration.x * TimeStep;
	d_ParticleArrayVelocity[gThID].y   += Acceleration.y * TimeStep;
	d_ParticleArrayVelocity[gThID].z   += Acceleration.z * TimeStep;

	// updating the coordinates
	d_ParticleArrayPosition[gThID].x += (d_ParticleArrayVelocity[gThID].x * TimeStep)/2.0;
	d_ParticleArrayPosition[gThID].y += (d_ParticleArrayVelocity[gThID].y * TimeStep)/2.0;
	d_ParticleArrayPosition[gThID].z += (d_ParticleArrayVelocity[gThID].z * TimeStep)/2.0;

	// Energy calculation
	d_EnergyArray[gThID].x = Acceleration.w; // PE
	d_EnergyArray[gThID].y = ( (d_ParticleArrayVelocity[gThID].x * d_ParticleArrayVelocity[gThID].x) +
	                           (d_ParticleArrayVelocity[gThID].y * d_ParticleArrayVelocity[gThID].y) +
							   (d_ParticleArrayVelocity[gThID].z * d_ParticleArrayVelocity[gThID].z) )
							 * d_ParticleArrayPosition[gThID].w;


}





//Auxiliary functions

void EnergyLog(char* FileName, REAL time, REAL Etot, REAL Ek, REAL Ep, REAL energyerr, char* fmode)
{

    //
    // used to generated the Energy Log.
    //

    FILE *file_p;

    file_p = fopen(FileName, fmode);

	fprintf(file_p, "%e %e %e %e %e\n", time, Etot, Ek, Ep, energyerr);

    fclose(file_p);

}


void LoadInitialConditionsData(char* FileName, unsigned int* N, REAL_CUDA** ParticleArrayPosition, REAL_CUDA** ParticleArrayVelocity, REAL2** EnergyArray)
{
    FILE* fptr;

    if ((fptr = fopen(FileName, "r")))
    {
        char    LineBuffer[BUFFER_MAX];
        unsigned int     LineCount = 0;

        // The first line of the initial conditions file should contain the amount of particles information.
        // So, read this value and use it for allocate memory for the particles arrays.
        fgets (LineBuffer, BUFFER_MAX, fptr);
        sscanf(LineBuffer, "%d ", &LineCount);
        *N = LineCount; // 'N' value finally discovered!

        // allocate memory for the Particle arrays.
        *ParticleArrayPosition = (REAL_CUDA*)malloc( sizeof(REAL_CUDA) * LineCount );
        *ParticleArrayVelocity = (REAL_CUDA*)malloc( sizeof(REAL_CUDA) * LineCount );
		*EnergyArray           = (REAL2*)malloc    ( sizeof(REAL2)     * LineCount );

        // read the file line by line, saving position and velocity into data arrays.
        // each file line represents only one particle.
        for (int i=0; i<LineCount; i++)
        {
            fgets (LineBuffer, BUFFER_MAX, fptr);
            sscanf(LineBuffer, "%lf %lf %lf %lf %lf %lf %lf " , &(*ParticleArrayPosition)[i].x,      // PosX ==> &(((*ParticleArrayPosition)+i)->x),
            //sscanf(LineBuffer, "%f %f %f %f %f %f %f " , &(*ParticleArrayPosition)[i].x,      // PosX ==> &(((*ParticleArrayPosition)+i)->x),
                                                         &(*ParticleArrayPosition)[i].y,      // PosY
                                                         &(*ParticleArrayPosition)[i].z,      // PosZ
                                                         &(*ParticleArrayVelocity)[i].x,      // VelX
                                                         &(*ParticleArrayVelocity)[i].y,      // VelY
                                                         &(*ParticleArrayVelocity)[i].z,      // VelZ
                                                         &(*ParticleArrayPosition)[i].w);     // Mass
			(*EnergyArray)[i].x = 0.0;
			(*EnergyArray)[i].y = 0.0;
			//((*EnergyArray)+i)->x = 0.0;
			//((*EnergyArray)+i)->y = 3.5;
        }

        //close the file and release ay used resource
        fclose(fptr);

        return;

    }
    else
    {
        //Error! Couldn't open the file!
        printf("\nError opening file %s!\n", FileName);
        return;
    }
}

void GenerateResultFile(char* FileName, REAL_CUDA* ParticleArrayPosition, REAL_CUDA* ParticleArrayVelocity, unsigned int N)
{

    //
    // used to generated the data file containing the updated particles positions.
    //

    FILE *file_p;

    file_p = fopen(FileName, "w");

    //first line must cointain the number of particles
    //fprintf(file_p, "%d\n", N);

    for(int i=0; i<N; i++)
    {
        //fprintf(file_p, "%f %f %f %f %f %f %f\n", ParticleArrayPosition[i].x, // teste com uso de double4
	    fprintf(file_p, "%e %e %e %e %e %e %e\n", ParticleArrayPosition[i].x,
                                                  ParticleArrayPosition[i].y,
                                                  ParticleArrayPosition[i].z,
                                                  ParticleArrayVelocity[i].x,
                                                  ParticleArrayVelocity[i].y,
                                                  ParticleArrayVelocity[i].z,
                                                  ParticleArrayPosition[i].w);

    }

    fclose(file_p);

}

void traceParticle(int id, REAL t, REAL* KE, REAL* PE, char* filemode)
{
	char ResultFileName[100];

	sprintf(ResultFileName, "./data/energy-particule%d.log", id);

	FILE *file_p;

    file_p = fopen(ResultFileName, filemode);

	fprintf(file_p, "%e %e %e\n", t, *KE, *PE);

    fclose(file_p);

    //soma das energias
    sprintf(ResultFileName, "./data/somaEnergy-particule%d.log", id);
    file_p = fopen(ResultFileName, filemode);

	fprintf(file_p, "%e\n",(REAL)(*KE+*PE));

    fclose(file_p);
}

int main(int argc, char **argv)
{
    printf("\nCUDASaganSimulator version 1.0 Beta\nStarting...\n:: Getting information about available GPU...\n");


	// Getting information about installed GPU
	// this information will be used to calculate the size of execution kernel blocks.
    int devID;
    int MAX_NUMBER_THREADS_PER_BLOCK;
    int NUMBER_OF_BLOCKS;
    int THREADS_PER_BLOCK;
    cudaDeviceProp props;
    CUDA_SAFE_CALL( cudaGetDevice(&devID) );
    CUDA_SAFE_CALL( cudaGetDeviceProperties(&props, devID) );
    if (props.major == 9999 && props.minor == 9999)
    {
        // there is not CUDA devices available!!!
        printf("\n ==>>> There is no CUDA device available! Exiting now...\n");
        return 0;
    }
    MAX_NUMBER_THREADS_PER_BLOCK = props.maxThreadsPerBlock;

    printf(":::: GPU found: %s\n", props.name);
    printf(":::: Maximum Number os threads per block allowed: %d\n", MAX_NUMBER_THREADS_PER_BLOCK);




    // discrete variables
    //char InitialConditionsFileName[50];
    char InitialConditionsFileName[] = "esfera395CUDA.ini";
    char strTrace[20];
    //int trace;
    char    Arquivo[100];
    unsigned int     N = 0;

    if (argc>=2)
    {
		int argcont;

		for(argcont=0;argcont<argc;argcont++)
		{

			if (!strcmp(argv[argcont],"-init"))
			{
				strcpy(InitialConditionsFileName,argv[argcont+1]);
			}

			if (!strcmp(argv[argcont],"-trace"))
			{
				strcpy(strTrace,argv[argcont+1]);
			}
		}
	}

	else
	{
		printf("\n:::Initial file not given, please start the program using -init filename.ini:::\n\n");
		exit(1);
	}

	sprintf(Arquivo,"./data/%s",InitialConditionsFileName);
	//if (strlen(strTrace)>0)
	//	trace = atoi(strTrace);
	//else
	//	trace = 0;

	// Auxiliary variables to monitor the energy of the n-body system.
	//REAL* traceKE;   // Kinetic Energy at host memory
	//REAL* tracePE;   // Potential Energy at host memory
	//REAL* d_traceKE; // Kinetic Energy at device memory
	//REAL* d_tracePE; // Potential Energy at device memory
	//traceKE = (REAL*)malloc(sizeof(REAL));
    //tracePE = (REAL*)malloc(sizeof(REAL));
    //cudaMalloc((void**)&d_traceKE, sizeof(REAL));
    //cudaMalloc((void**)&d_tracePE, sizeof(REAL));

    //REAL flag = 0;
    //REAL* d_flag;


    // host memory for particle arrays.
    REAL_CUDA* ParticleArrayPosition;
    REAL_CUDA* ParticleArrayVelocity;

	//host memory for energy monitoring.
	REAL2* EnergyArray;

    // Load 'N' value and particle information from initial conditions file.
    printf(":: Loading Initial Conditions File (%s)...\n", Arquivo);
    LoadInitialConditionsData(Arquivo, &N, &ParticleArrayPosition, &ParticleArrayVelocity, &EnergyArray);

    // Once 'N' is now knowed, it's possible to define the size of the execution block
    printf(":: Initial Conditions File loaded. %d particles were found.\n", N);

	/*for(int u =0; u<N; u++)
	  printf(" EnergyArray[%d].x = %f; EnergyArray[%d].y = %f\n", u, EnergyArray[u].x, u, EnergyArray[u].y);

	return 0;
	*/

    /*
    if(N>MAX_NUMBER_THREADS_PER_BLOCK)
    {
        NUMBER_OF_BLOCKS =  ((N % MAX_NUMBER_THREADS_PER_BLOCK)>0 ? 1 : 0) + (N / MAX_NUMBER_THREADS_PER_BLOCK);
        THREADS_PER_BLOCK = MAX_NUMBER_THREADS_PER_BLOCK;
    }
    else
    {
        NUMBER_OF_BLOCKS = 1;
        THREADS_PER_BLOCK = N;
    }
    */
	float cc;
	cc = props.major+props.minor*0.1;
   	printf("\nCC: %f",cc);

	int THREADS_KERNEL;

	if (cc <1.3)
		THREADS_KERNEL = 16;
	else if (cc==1.3)
			THREADS_KERNEL = 26;
	else
			THREADS_KERNEL = 34;

	printf("\nThreads kernel: %d",THREADS_KERNEL);

	int warps_block = floor(props.regsPerBlock/(THREADS_KERNEL*props.warpSize));
	THREADS_PER_BLOCK = warps_block * props.warpSize;

	if (THREADS_PER_BLOCK>MAX_NUMBER_THREADS_PER_BLOCK)
	{
		THREADS_PER_BLOCK = MAX_NUMBER_THREADS_PER_BLOCK;
	}

	//este vai ser o numero maximo sem estourar o numero de regs
	MAX_NUMBER_THREADS_PER_BLOCK = THREADS_PER_BLOCK;

	int it;
	int warps_per_sm;
	int warps_per_block;
	int regs_per_block;
	int blocks_per_sm_regs;
	int blocks_per_sm_warps;
	int blocks_per_sm_shared;
	int warps_alloc_granularity;
	int rounded_warps_per_block;
	int regs_alloc_unit_size;
	int max_warps_per_sm;
	int max_blocks_per_sm;

	//set warps_alloc_granularity, values from nvidia
	if (cc >= 3.0)
	{
		warps_alloc_granularity = 4;
		max_blocks_per_sm = 16;
	}
	else
	{
		warps_alloc_granularity = 2;
		max_blocks_per_sm = 8;
	}

	//set regs_alloc_unit_size and max_warps_per_sm, values from nvidia
	if (cc < 1.2)
	{
		regs_alloc_unit_size = 256;
		max_warps_per_sm = 24;
	}
	else if (cc < 2)
	{
		regs_alloc_unit_size = 512;
		max_warps_per_sm = 32;
	}
	else if (cc< 3)
	{
		regs_alloc_unit_size = 64;
		max_warps_per_sm = 48;
	}
	else
	{
		regs_alloc_unit_size = 256;
		max_warps_per_sm = 64;
	}



	//max_threads é o numero de threads por bloco com maior ocupação
	//max_warps_per_mp é o numero de warps que o max_threads vai ocupar do mp
	int best_warps_per_sm=0;
	int best_threads=0;

	/*ARRUMAR REGS_PER_BLOCK */
	for(it=32;it<MAX_NUMBER_THREADS_PER_BLOCK;it+=32)
	{
		warps_per_block = it/props.warpSize;

		if (cc<2.0) //so allocation granularity = block
		{
			printf("\nallocation granularity = block");
			if (warps_per_block%warps_alloc_granularity != 0) //arredondar para cima warps_per_block
			{
				rounded_warps_per_block=(warps_per_block/warps_alloc_granularity)*warps_alloc_granularity+warps_alloc_granularity;
				regs_per_block=rounded_warps_per_block*THREADS_KERNEL*props.warpSize;
				if (regs_per_block%regs_alloc_unit_size != 0) //arredondar para cima regs_per_block
					regs_per_block = (regs_per_block/regs_alloc_unit_size)*regs_alloc_unit_size+regs_alloc_unit_size;
			}
		}
		else //so allocation granularity = warps
		{
			printf("\nallocation granularity = warps");
			if (THREADS_KERNEL*props.warpSize%regs_alloc_unit_size != 0) //arredondar para cima THREADS_KERNEL*props.warpSize
			{
				if (warps_per_block%warps_alloc_granularity != 0) //arredondar para cima warps_per_block
					regs_per_block = ((THREADS_KERNEL*props.warpSize/regs_alloc_unit_size)*regs_alloc_unit_size+regs_alloc_unit_size)*((warps_per_block/warps_alloc_granularity)*warps_alloc_granularity+warps_alloc_granularity);
				else
					regs_per_block = ((THREADS_KERNEL*props.warpSize/regs_alloc_unit_size)*regs_alloc_unit_size+regs_alloc_unit_size)*warps_per_block;
			}
			else //nao precisei arrendondar THREADS_KERNEL*props.warpSize
			{
				//OK
				if (warps_per_block%warps_alloc_granularity != 0) //arredondar para cima warps_per_block
					regs_per_block = THREADS_KERNEL*props.warpSize*((warps_per_block/warps_alloc_granularity)*warps_alloc_granularity+warps_alloc_granularity);
				else //nao precisei arrendondar warps_per_block
					regs_per_block = THREADS_KERNEL*props.warpSize*warps_per_block;
			}
		}

		//blocos por sm devido ao numero de regs
		blocks_per_sm_regs = props.regsPerBlock/regs_per_block;
		//printf("\nprops.regsPerBlock: %d, regs_per_block: %d, blocks_per_sm_regs: %d",props.regsPerBlock,regs_per_block,blocks_per_sm_regs);

		//blocos por sm devido ao numero de warps
		blocks_per_sm_warps = max_warps_per_sm/warps_per_block;
		if (blocks_per_sm_warps>max_blocks_per_sm) //minimo deles
			blocks_per_sm_warps=max_blocks_per_sm;

		//blocos por sm devido ao numero de shared memory
		blocks_per_sm_shared = props.sharedMemPerBlock/SHARED_MEMORY_SIZE;

		if (blocks_per_sm_regs <= blocks_per_sm_warps && blocks_per_sm_regs <= blocks_per_sm_shared)
			warps_per_sm = blocks_per_sm_regs * warps_per_block;
		else if (blocks_per_sm_warps <= blocks_per_sm_shared)
			warps_per_sm = blocks_per_sm_warps * warps_per_block;
		else
			warps_per_sm = blocks_per_sm_shared * warps_per_block;

		//printf("\nblocks_per_sm_regs: %d, blocks_per_sm_warps: %d, blocks_per_sm_shared: %d, warps_per_sm: %d",blocks_per_sm_regs,blocks_per_sm_warps,blocks_per_sm_shared,warps_per_sm);

		if (warps_per_sm > best_warps_per_sm && (N/it + N%it>0?1:0)<=props.multiProcessorCount*max_blocks_per_sm)
		{
			best_warps_per_sm = warps_per_sm;
			best_threads = it;
		}
	}

	THREADS_PER_BLOCK = best_threads;

	//teste
	//THREADS_PER_BLOCK = 512;

	if (N <= THREADS_PER_BLOCK)
	{
		NUMBER_OF_BLOCKS = 1;
	}

	else
	{
		NUMBER_OF_BLOCKS =  (N % THREADS_PER_BLOCK > 0 ? 1 : 0) + N / THREADS_PER_BLOCK;
	}


    printf(":: Number of execution blocks = %d\n", NUMBER_OF_BLOCKS);
    printf(":: Number of threads per block = %d\n", THREADS_PER_BLOCK);



    // device memory allocation for particle data arrays.
    CUDA_SAFE_CALL ( cudaMalloc((void**)&d_ParticleArrayPosition, sizeof(REAL_CUDA)*N) );
    if ( cudaSuccess != cudaGetLastError() )
	    printf( "Error malloc position!\n" );

    CUDA_SAFE_CALL ( cudaMalloc((void**)&d_ParticleArrayVelocity, sizeof(REAL_CUDA)*N) );
    if ( cudaSuccess != cudaGetLastError() )
	    printf( "Error malloc velocity!\n" );

    if (d_ParticleArrayPosition ==NULL || d_ParticleArrayVelocity == NULL)
    {
		printf("\nErro no cudaMalloc!");
		exit(1);
	}

	// device memory allocation for energy monitoring array.
    CUDA_SAFE_CALL ( cudaMalloc((void**)&d_EnergyArray, sizeof(REAL2)*N) );
    if ( cudaSuccess != cudaGetLastError() )
	    printf( "Error malloc Energy Array on device memory!\n" );

    if (d_EnergyArray == NULL)
    {
		printf("\nErro no cudaMalloc do EnergyArray!\n");
		exit(1);
	}


    // copy particle information already loaded from initial conditions file to data arrays on device memory.
    CUDA_SAFE_CALL(cudaMemcpy(d_ParticleArrayPosition, ParticleArrayPosition, sizeof(REAL_CUDA)*N, cudaMemcpyHostToDevice));
    if ( cudaSuccess != cudaGetLastError() )
	{
		printf("erro no cudamemcpy");
	}

    CUDA_SAFE_CALL( cudaMemcpy(d_ParticleArrayVelocity, ParticleArrayVelocity, sizeof(REAL_CUDA)*N, cudaMemcpyHostToDevice) );
	if ( cudaSuccess != cudaGetLastError() )
	{
		printf("erro no cudamemcpy");
	}

    CUDA_SAFE_CALL( cudaMemcpy(d_EnergyArray, EnergyArray, sizeof(REAL2)*N, cudaMemcpyHostToDevice) );
	if ( cudaSuccess != cudaGetLastError() )
	{
		printf("erro no cudamemcpy");
	}


    int i=0;    // we need these sequential numbers to name the output files.
    REAL t = 0.0;        // controls the amount of time steps already processed for the simulation.
    REAL dt = 0.0001;    // indicates the contribution of each scan for the simulation evolution.
    REAL dtout = 0.01;   // indicates the value of the step where a result file has to be generated.
    REAL tout = t+dtout;
    int    timeSteps = 10; //WARNING!! Don't forget to update this value with the desired time steps amount!
    REAL epson = 0.005;
    if(N>100)epson = 4.0/N; // softening factor


    char ResultFileName[100];


	//
	// Definindo o tamanho da sharedmemory a ser alocada. queremos limitar em 16Kb por bloco
	//
	int tileWidth = 256; // quantidade de posiçoes a serem alocadas na shared memory por bloco.
	// assumindo sizeof(double) = 16bytes, logo double4 = 16 * 4 = 64bytes
	// entao, 256 posicoes de double4 = 16K
	size_t SharedMemorySize = tileWidth * sizeof(REAL_CUDA);

	// help to monitor the energy of the particle system
	bool flagCtrlEnergyError = true;
	REAL InitialTotalEnergy = 0.0;	

	// time counting
	cudaEvent_t Start;
	cudaEvent_t Stop;
	float ExecutionTime;
	cudaEventCreate(&Start);
	cudaEventCreate(&Stop);

	cudaEventRecord(Start, 0);

	printf("\n::: processing...\n");
    while (t<=timeSteps)
    {

		// run the main kernel
		MainKernel<<<NUMBER_OF_BLOCKS, THREADS_PER_BLOCK,SharedMemorySize>>>(d_ParticleArrayPosition, d_ParticleArrayVelocity, d_EnergyArray, N, epson, dt, tileWidth);


       cudaError err = cudaGetLastError();
        if ( cudaSuccess != err )
        {
			printf("\n!!! Error running MainKernel !!!\n ::: File: %s\n ::: Line: %d\n ::: CUDA Error: %s\n", __FILE__, __LINE__, cudaGetErrorString(err));
			//getchar();
		}


        t += dt;

        if(t>tout)
        {
                // copy the updated data arrays from device global memory back to host memory.
                CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayPosition, d_ParticleArrayPosition, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );
                CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayVelocity, d_ParticleArrayVelocity, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );

                // copy auxiliary energy monitoring variables back for host memory.
 			    CUDA_SAFE_CALL( cudaMemcpy(EnergyArray,   d_EnergyArray,   sizeof(REAL2)*N, cudaMemcpyDeviceToHost) );

                tout = t + dtout;
                sprintf(ResultFileName, "./data/%03d.dat", i++);

                //printf("\n::: generating result file %s...\n", ResultFileName);
                if(i%10==0)
                {
					printf(" __%d__ ", i);
				}
				else
				{
					printf(".");
				}
                GenerateResultFile(ResultFileName, ParticleArrayPosition, ParticleArrayVelocity, N);

                // Energy Logging
                REAL EK = 0.0;
                REAL EP = 0.0;
				for(unsigned int i=0; i<N; i++)
				{
					EP += EnergyArray[i].x;
					EK += EnergyArray[i].y;
				}

                REAL TotalEnergy  = EK/2.0 + EP;
                
                REAL EnergyError = 0.0; 
                
                if(flagCtrlEnergyError) 
                {
					// if it's the first iteraction, takes the energy returned by kernel as initial energy for the particle system.
					InitialTotalEnergy = TotalEnergy;
					flagCtrlEnergyError = false;
					
					EnergyLog("./data/energy.log", t, TotalEnergy , EK, EP, 0.0, "w");
					//printf("\n ::: InitialTotalEnergy = %e\n", InitialTotalEnergy);
					//getchar();
				} 
				else 
				{
					EnergyError = ( TotalEnergy - InitialTotalEnergy) / fabs(InitialTotalEnergy);                
				
					EnergyLog("./data/energy.log", t, TotalEnergy , EK, EP, EnergyError, "a");
					//printf("\n ::: EnergyError = %e :: TotalEnergy = %e :: InitialTotalEnergy = %e :: fabs(InitialTotalEnergy) = %e \n", EnergyError, TotalEnergy, InitialTotalEnergy, fabs(InitialTotalEnergy) );
				}


        }

    }
    
    cudaEventRecord(Stop, 0);
    cudaEventSynchronize(Stop);
    
    cudaEventElapsedTime(&ExecutionTime, Start, Stop);
    
    printf("\n:::: Execution time: %f ms\n", ExecutionTime);


    // copy the updated data arrays from device global memory back to host memory.
    CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayPosition, d_ParticleArrayPosition, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );
    CUDA_SAFE_CALL( cudaMemcpy(ParticleArrayVelocity, d_ParticleArrayVelocity, sizeof(REAL_CUDA)*N, cudaMemcpyDeviceToHost) );

    cudaFree(d_ParticleArrayPosition);
    cudaFree(d_ParticleArrayVelocity);
    free(ParticleArrayPosition);
    free(ParticleArrayVelocity);



    printf("\nValor de N: %d\n", N);


    printf("\n\nYou ran me by typing: %s\n", argv[0]);
    printf("\nPlaca: %s\n", props.name);
    return (0);
}

